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Abstract. We study transport through an interacting model system consisting of a 
central correlated site coupled to finite bandwidth tight-binding leads, which are considered 
as effectively noninteracting. Its nonequilibrium properties are determined by real-time 
propagation of the Kadanoff-Baym equations after applying a bias voltage to the system. 
The electronic interactions on the central site are incorporated by means of self-energy 
approximations at Hartree-Fock, second Born and GW level. We investigate the conditions 
under which multiple steady-state solutions occur within different self-energy approximations, 
and analyze in detail the nature of these states from an analysis of their spectral functions. 
At the Hartree-Fock level at least two stable steady-state solutions with different densities and 
currents can be found. By applying a gate voltage-pulse at a given time we are able to switch 
between these solutions. With the same parameters we find only one steady-state solution when 
the self-consistent second Born and GW approximations are considered. We therefore conclude 
that treatment of many-body interactions beyond mean-field can destroy bistability and lead to 
qualitatively different results as compared those at mean-field level. 



1. Introduction 

The experimental observation [HE] of a hysteresis loop in the I/V characteristic of double-barrier 
resonant tunneling structures prompted intense theoretical activities to gain a microscopic 
understanding of this phenomenon. Several authors have been able to reproduce the hysteresis 
behavior by treating the Coulomb interaction at a mean field level (3J IU EJ |6]. Self- 
consistent calculations have revealed the presence of bistable solutions, one of the solution 
being characterized by a considerable accumulation of charge in the potential well. Subsequent 
experimental work on double-barrier resonant tunneling diodes has, however, demonstrated that 
hysteresis loops do not always occur [7|. The suppression of the intrinsic bistability has been 
attributed to exchange-correlation effects jS]. 



With the advent of molecular electronics [§] the study of intrinsic bistability in nanoscale 
devices has gained attention due to the possibility of developing, for example, molecular diodes. 
So far, most of the work has focused on the steady-state I/V curve of molecular junctions 
attached to metallic leads. Up to date calculations are performed within the one-particle 
scheme of time- independent density functional theory (DFT). At the Hartree level bistability 
was reported by Negre el. al for a double quantum dot structure [TO] . 

The mechanism of bistability and the calculation of switching times between two different 
states are mostly unexplored and the question how correlations affect the bistability has received 
very scarce attention. The purpose of the present exploratory paper is twofold: to extend the 
analysis to the time-domain and to study the role of memory effects in a bistable interacting 
resonant level model (IRLM). 

Two complementary theoretical approaches will be used for calculating the time-dependent 
current and density, namely Time-Dependent (TD) DFT and Many-Body Perturbation Theory 
(MBPT). TDDFT [11] provide an exact framework to account for correlation effects both in the 
leads and the central region [12] . Within TDDFT the basic quantities that are propagated in time 
are the one-particle orbitals which depend on only one time variable. This property renders the 
implementation computationally favorable [13] . Most approximations to the TDDFT potential, 
however, do not include memory effects and the exchange-correlation part is approximated by 
local or semi-local functionals of the density. The lack of more sophisticated approximations 
represents, at present, a major obstacle for an accurate first principle description of TD 
quantum transport through interacting regions. MBPT has the advantage over TDDFT of 
allowing for a systematic inclusion of relevant physical processes through a selection of Feynman 
diagrams. Thus, MBPT provides an important tool to proceed beyond the commonly used 
adiabatic approximations and to quantify the importance of memory effects through advanced 
approximations to the self-energy. We recently proposed [14} [T5] a time dependent MBPT 
formulation of quantum transport, based on the real-time propagation of the Green function 
|16 [ ll7 l fT8] for open and interacting systems. First the Dyson equation for the connected system 
is solved to self-consistency on the imaginary axis. After that the Green function is propagated 
with the Kadanoff-Baym equations using different level of conserving approximations. As the 
Green function depends on two time variables the implementation of the MBPT scheme is more 
demanding than that of the TDDFT scheme. We expect that the interplay between MBPT and 
TDDFT will be essential to develop accurate approximation for systems more complex than the 
one studied here. 

2. Interacting Resonant Level Model 

We study an Anderson-type of system [19] where the impurity is an interacting site coupled to 
the infinite one dimensional tight-binding leads of finite band width. The Hamiltonian of the 
system reads 

(J a,a' i.a,cr 

i,a,cr a, it 

where i denotes the site indices and a is the spin index, eq is the on-site energy of the localized 
site, IA is the strength of the two-particle interaction on the central site, b is the hopping 
parameter between lead sites, U a (t) is time-dependent bias voltage in the leads (a = L/R), 
a is the on-site parameter in the leads and Vo t i a denotes coupling between the lead and the 
localized site. The fermionic creation- and annihilation operators in the leads a are denoted 



as & , c whereas for the localized site they are denoted as aft , d. The quantity V g (t) denotes a 
time-dependent gate voltage. 

The main reason for studying the IRLM is that for this system the multiple steady-state 
solutions are easily found from a fixed-point equation for the density on the localized site. 
The IRLM is the simplest system in which bistability occurs and hence allows for a clear 
interpretation of its multiple steady-state solutions. We study the system at a finite temperature 
and under a finite bias, such that we are out of the Kondo regime in which the IRLM is often 
used. 
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Figure 1: a) A schematic representation of the studied system. The correlated central region (C) 
is coupled to infinite one dimensional left (L) and right (R) tight-binding leads via a coupling 
Hamiltonian. b) The Keldysh contour c. Times on the lower branch t + are later than times on 
the upper branch t_. The imaginary track extends up to the inverse temperature f3. 



3. Theoretical approach 

3.1. Kadanoff-Baym equations 

We study the nonequilibrium properties of the IRLM by means of time-propagation of the 
Kadanoff-Baym equations for the nonequilibrium Green function [20]. We assume the system to 
be contacted and in equilibrium at a chemical potential // and at inverse temperature (3 before 
the time t = Iq. For t > to the system is driven out of equilibrium by an external bias and 
we aim to study the time-evolution of the electron density and current. We give here a brief 
description of the approach described in more detail in Ref. |14j . 

The Keldysh Green function is defined as the expectation value of the contour ordered 
product [17] 

m 2) = , TV{L>(to-i/3,to)r c [^(l)^(2)]} 
Tr{l>(t - i/Mo)} 



where 7c denotes the time ordering operator along the Keldysh contour c (see Fig. lb ) and U is 
the time evolution operator. The average is taken over the grand canonical ensemble. We use the 
compact notation 1 = (xi,ii) and 2 = (x2,i2) where x = (r,c) is a collective space-spin index. 
From the Green function it is possible to calculate any one particle property of the system. For 
example the time-dependent density is given as 

<n(l)) = -zG<(l,l+), (3) 



where t + approaches t from an infinitesimally later time t + = t + 5. The equation of motion for 
the full system can be derived from the definition of the Green function Eq. Q and reads 



[id h -if(l)]G(l,2) = 6(1,2) + J d3E MB [G](l,3)G(3,2), (4) 

where H is the one-body part of the Hamiltonian. The self-energy £[G] incorporates all of 
the effects of exchange and correlation in the central region and is a functional of the Green 
function |17} IT6] . In a localized basis the one-body part of the Hamiltonian ([I]) and the many- 
body self-energy can be written in a block-matrix form 

( H LL H LC \ /0 0\ 

H=\Hcl H CC H CR , E MB = £$£[G C c] , (5) 
\ H RC H rr ) \0 0/ 

where the H aa and Hqc components describe the leads and the central system respectively, 
whereas the off-diagonal components describe the hopping between them |14j . We only 
consider the central region as interacting whereas the leads are effictively noninteracting. As 
a consequence, the many-body self-energy in Eq. ^ has non-vanishing elements only for the 
central region because the diagrammatic expansion starts and ends with an interaction line. In 
this work the electronic interactions are incorporated in £ MB [G] at HF, 2B and GW level [H]. 

Solving the problem of an open infinite system is equivalent to solving the problem of a closed 
system with an equation of motion which considers the leads through an embedding term |14j . 
In our case this reads (in the reminder of this paper we will suppress the spatial indices of the 
objects involved) 

[idt - H cc (t)] G CC (t, = S(t, t') +Jdi{ [E em (i, t) + X%g[Gcc]{t, t)] G CC (t, t')} , (6) 

where the embedding self-energy T, em (t, t') accounts for the tunneling of electrons between leads 
and central region. In its general form, the embedding self-energy reads [14| fT5| I2T] 

£em(M0 = y^em.aCMO = ^ Hca9aa(t,t')H a c, (7) 
a a 

where g a a(t,t') is the lead Green function and Hc a ,H a c is the coupling Hamiltonian between 
the central site and the leads. 

The current through the lead a can be expressed in terms of Keldysh Green functions as 

nuns] 

l a (t) = 2Re{ Tr c [ I dt[G< c (t, t)E^ a (t, t) + I diG% c (t, t)S< mjQ (t, t)} 

dfG 1 co (t,f)EL, a (f,t)l), (8) 



where we integrated on the Keldysh contour and where the superscripts A, R, < refer to 
advanced/retarded/lesser component of Green function/self-energy and ],[ are the mixed 
components having one time argument on a imaginary axis and the other on the real axis 
|20| \TE[ I14j . The trace is taken over the central region. The current accounts for the initial 
many-body and embedding effects through the last term in equation (TsT) which is an integral 



over the vertical track. Equation ([8]) is a generalization of the Meir-Wingreen formula |22j . We 
further define the nonequilibrium spectral function 

A(7» = -ImTr| ^e™[G> - G<](T + T -, T - T -\ (9) 

where r = t—t' is a relative time and T = (t+t')/2 is an average time coordinate. In equilibrium, 
this function is independent of T and has peaks below the Fermi level at the electron removal 
energies of the system, while above the Fermi level it has peaks at the electron addition energies. 
If the time-dependent external field becomes constant after some switching time, then also the 
spectral function becomes independent of T after some transient period and has peaks at the 
addition and removal energies of the biased system. 



3.2. Time- dependent density functional theory 

TDDFT provides an alternative framework to describe electron transport through interacting 
systems. In TDDFT [23] the time-dependent density of the interacting system is obtained 
through time-propagation of a Kohn-Sham system in an effective local potential. While in 
MBPT the correlation level is determined by the choice of the many-body self-energy, in TDDFT 
the main approximation is the functional used for the effective potential. The difficulty of 
TDDFT, compared to MBPT is the current lack of sufficient accurate approximations to the 
time-dependent exchange-correlation potential. However, the computational effort is much lower 
compared to a MBPT propagation. 

The problem brought forward by considering an open system can be resolved in a very similar 
manner as in MBPT, with the help of an embedding self-energy. The equation of motion for 
the fc-th single-particle orbital projected onto the central region, tl>k,c(t), reads 

[idt-Ecc(t)]1> k , c (t) = f dfS^Ct^Vfc.c^ + J^^^C*' )^^). ( 10 ) 

Jo 

where T,^ m (t,t) is the retarded embedding self-energy (see Eq. Q) and g^ a is the retarded lead 
Green function. The time-dependent density in the central region is obtained by 

occ 

n(t)=J2\Mt)\ 2 > (11) 

k 

where the summation is taken over all occupied orbitals in the time-dependent Slater determinant 



[24] . The technique to propagate Eq. (10) is described in detail in Ref. [13 . In this work we 
used this approach at a level, in which, for the system studied, the local exchange potential is 
equal to half the Hartree potential of the localized site. The results were found, as expected, to 
be identical to those obtained from the Kadanoff-Baym approach at HF level. 



3.3. Steady-state density 

We begin our study of the bistable regime by deriving an equation for the density on the localized 
site. This quantity is given by the lesser Green function at equal times 

(h(t)) = -iG<(t,t + ). (12) 



Since we consider the steady-state, we can assume that in the long-time limit the Green functions 
depend only on the relative time coordinate t — t'. In that case the Green function can be Fourier 



transformed with respect to the relative time variable and the expression for the steady-state 
becomes 

.-/&o<M. (13) 

The Green function appearing in this expression satisfies the equation 

G<{u) = G R {u)?.< t {u)G A {u), 



(14) 



where T,f ot {ui) = Y,f m {uj) + Y,^ R)< (uj) and where G A = [G R ] . The retarded Green function has 
the expression 

(15) 



G R {u) = [u-e - Re[£* ( w )] - iIm[S« (w)]] 1 



For the tight-binding leads, that we consider, the retarded embedding self-energy is given by 
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(16) 



where ui a = uj — a + [i — U a , with the lead-on-site parameter a and the hopping parameter 
between the lead sites b. The chemical potential is denoted by /J,, the applied bias for the lead a 
by U a , and Vi a> o, Vo,1q are the left/right couplings between leads and the central site. The lesser 
component of the embedding self-energy can be expressed as T,f m<Q ,(uj) = if a (u)T a (uj) where f a 



is the Fermi distribution of lead a and r Q is defined in Eq. (16) 



If we integrate the left hand side of Eq. (14) over all frequencies then according to Eq. (13) 



we obtain an expression for the density per spin n on the localized site 
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where A = A^j + Al and T = Tr + Tl. This is a Meir-Wingreen-type equation for the 
density |25} [22] and is valid for transport through interacting systems. Within the HF 
approximation for the IRLM one has T, MB ^< = Im[S MB ' K ] = and Re[T, MB ' R ] = Un. Note that 
we include the time-singular part of the self-energy in the definition of the retarded/ advanced 
component, see e.g., [18]. Within this approximation, Eq. (17) now becomes a self-consistent 



fixed-point equation for the density n, since the value of the integral ( 17) depends on the density 
via the term U n in the denominator. 



4. Results and discussion 

We consider a biased system with the following parameters: Vo,lr = Vo,il = V = —0.35, Ul 
= 1.5, U R = 0.0, U = 2.0, b = -0.5, a = // and fx = 0.3, = 90. The leads are half-filled 
such that the Fermi level of lead a is positioned at \x + U a and the width of the lead band is 
[fi + U a — 26, \x + U a + 26]. With these parameters, within the HF approximation, the Meir- 
Wingreen approach [251 E2] predicts three solutions for the steady-state density n in Eq. (17). 
The three fixed points are shown in the inset of Fig. 2a, where we display the left and right hand 
side of Eq. (17). The corresponding densities are n\ = 0.33, n<i = 0.58 and 713 = 0.66, which 
should be compared to the density of no = 0.28 of the unbiased equilibrium state. 

In 2a we show the steady-state spectral functions corresponding to these three solutions, 
which can be obtained directly from the retarded Green function G r {oj). The peak of the 
spectral function corresponding to density n\ = 0.33 is positioned at the energy 0.5. We thus 
see that this spectral function is located in an energy range within the right lead energy-band 
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Figure 2: a) Spectral functions for the HF approximation corresponding to different steady-state 



solutions. Inset: The graphical solution of the integral in Eq. (17): (1) is the left hand side of 
the equation whereas (2) is the value of the integral on the right hand side with given density 
n. b) Steady-state spectral functions for the 2nd Born and GW approximations. 



(see right side of Fig. 2a). The spectral function corresponding to density 77,3 = 0.66 is peaked 
approximately at energy 1.5 and is located in an energy range within the left lead energy-band. 
The HF spectral function corresponding second solution with density (77,2 = 0.58) is peaked 
on the top edge of the right lead energy-band and located in between the spectral functions 
corresponding to the densities n\ and n%. By time-propagation (see below), we find that the two 
solutions corresponding to densities n\ and 77,3 lead to stable steady-states, i.e, states that are 
reachable by time-propagation after applying a bias. On the other hand, the state corresponding 
to density 712 is unstable and cannot be reached by time-propagation. The spectral function 
corresponding to the state with density 77,2 has a large overlap with the one of density 773. This 
indicates that (for sufficiently slow gate switching, see below) during the time-propagation more 
charge will flow onto the central site resulting in a stable steady-state with density 71.3. From 
the analysis of the spectral functions we thus conclude that the density bistability in this system 
occurs when the spectral functions of the two stable solutions are localized well within one of 
the lead energy-bands and are well-separated. This happens exactly when the leads have a small 
overlap and the system is within the region of negative differential resistance (NDR). 



Let us now consider the situation when we go beyond the HF approximation. In Fig. 2b 
we show the steady-state spectral functions for the 2B and GW approximations, as obtained 
from time-propagation of the Kadanoff-Baym equations. With these approximations we found 
only one steady-state solution, with very broad spectral function due to enhanced quasi-particle 
scattering at finite bias |14j . The spectral weight of the 2B and GW states is spread almost 
uniformly over the energy range from the bottom of the right lead energy band to the top of left 
lead energy band and extends well outside the lead bands. We also observe two small peaks in 
these spectral functions which occur approximately at 0.8 and 1.3, corresponding to the edges 
of the lead energy bands. 

We now go to the time domain and consider how the steady-states at HF level, corresponding 
to densities ri\ = 0.33 and 713 = 0.66, can be reached by time-propagation starting from the 
initially unbiased equilibrium state. We generate these steady-states by applying time-dependent 
gate pulses and we use the same technique to switch between them. In this work we have used 



an exponentially decaying gate voltage of the form 



Vg(t) 



r V a e* 

_ Vge -l(t-T g ) 



if < t < T g 

if Tg < t < 2Tg 

if t > 2T n 



(18) 



The steady-state with density n\ is obtained by time-propagation after applying a sudden 
constant bias f7z,(i) = TJzfi{£) in the left lead, without applying the gate voltage. This is displayed 



in Fig. 3a The steady-state of highest density 713 is obtained (in addition to switching on the 



sudden bias in the leads) by applying an exponentially decaying gate voltage to the impurity 
site (with amplitude V g = 1.5, decay rate 7 = 0.1 and T g = 00). 
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Figure 3: a) Two different steady-state densities and the switching between the solutions by 



applying a gate voltage in the form of the equation (18). b) Two different steady-state currents 
through right lead and the switch between the solutions with a gate voltage of the form (18). 



In Fig. 3a we show the time evolution of the density for various switchings between the 



steady-states. After we apply the first (sufficiently slow) gate voltage of equation (18) we let 
the system reach the steady-state of density 713. At a time T g we apply a second gate voltage in 
the opposite direction. The system shows a transient behavior after which, the system reaches 
the steady-state of density n\. If we apply the third gate voltage at time 2T g (see equation 
( |18| ) ) , the system exhibits a short transient behavior and attains again the original steady-state 
of density 723 . 

Corresponding to the densities m and 713 there are two distinguishable solutions for the 
currents In(t) flowing into the right lead. These are shown in Fig. 3b where the lower value of 
the current corresponds to the state of density n\ and the higher value of the current corresponds 
to the state of density 773 . The various transients observed in the density correspond directly to 
the transients in current shown in Fig. 3b The frequency of the strongly damped oscillations 



of the transients (observed upon switching between the states in Fig 3b ) is approximately given 
by the gate voltage, which causes a temporary change of V g in the level position. As one can see 
from the Fig. 3a if the decay rate 7 > 0.5, i.e., when the switching is too fast, the steady-state of 
density 773 = 0.66 cannot be reached from the initial ground state of density no = 0.28, because 
the system does not have enough time to acquire sufficient density. 

In Fig. 4o we show, within the HF approximation, the nonequilibrium spectral function 
A(T,lo) of Eq.Q for a switch (with T g = 60, V g = 1.5, and 7 = 0.1) from the steady-state with 



density n\ to the state with density 723. When the gate is applied, the upper side of the spectral 
peak, starts to oscillate. The spectral peak undergoes a transition and overshoots before it 
settles to a new value of 1.5 within the right lead band. During this transition we can observe 
the appearance of a sharp peak in the spectral function located approximately at energy 1.3 
corresponding to the top edge of the right lead band. This is shown in the inset of Fig. |4a[ 
where we display a snapshot of A(T, uj) at time T = 70 at which the central site has density 
n(T) = 0.53, which is very close to the density of unstable state U2- Therefore, although this 
state does not lead to a steady-state, it can still be observed in the nonequilibrium spectral 
function. 
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Figure 4: a) Nonequilibrium spectral function A(T, oj) within the HF-approximation for the 
switch from the density m to density n^. Inset: Snapshot of spectral function corresponding to 
density n%. b) The densities and currents calculated within the 2B and the GW approximations 
with and without an exponentially decaying gate voltage. 



In Fig. 4b we show the densities and the currents obtained within the 2B and the GW 



self-energy approximations. For these cases only one steady-state is obtained with a density 
of about 0.5 on the central site. We have applied a time-dependent gate voltage of the form 
V g (t) = Vge~ l1 for t > 0. If no gate is applied the transient regime is shorter and the density 
attains its steady-state value faster, compared to the case where the gate is applied. One can 
observe that the steady-state values of the currents and the densities for both GW and 2B are 
close to each other implying that for this system the single-bubble diagram, common to both 
approximations, plays a crucial role |14| . 

We conclude that for both the 2B and the GW approximation the spectral functions are very 
broad (see Fig. 2b ) which makes it impossible to locate these spectral functions within an energy 
range in one of the lead energy-bands leading to two well-separated states. As a consequence the 
bistable regime is lost in the 2B and GW approximations, at least for the parameters considered. 



5. Conclusions 

In this paper we have investigated the switching between the steady-states of an interacting 
resonant level model connected to leads. For a given set of parameters we used a Meir- 
Wingreen type of approach to identify three different steady-state values for the density within 
the HF approximation. We showed that by applying an exponentially decaying gate voltage 
pulse during the time-propagation, at Hartree-Fock level, we can reach two different solutions 
and switch between them. However, due to a strong quasi-particle broadening of the spectral 



function, bistability is lost when the second Born and GW approximations are considered. 
Hence treatment of many-body interactions beyond mean-field can destroy bistability and lead 
to qualitatively different results as compared to those at mean-field level. 
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